Pyromania—Life after a fiery death: fire and plant biomass loading affect dissolved organic matter in experimental ponds
CJ Spiegel et al. 2023, accepted at Global Change Biology
11/06/2023
Project overview
The goal of the Pyromania project is to test how terrestrial subsides (plant biomass loading) and burning influence DOM composition and degradation. We used a manipulative experiment to assess a range of plant material quantities (0-400g per tank) and fire treatment (burned vs unburned material) and the non-linearity of these effects on aquatic systems through 4 time-point samplings. We used 400L aquatic mesocosms and ran the experiment for ~90d in 2021-2022.
DATA SETS
This data set is among 3 to be generated for the project and focuses
on:
- dissolved organic carbon (DOC) and
- Fluorescence and UV-vis spectroscopic absorbance
indices
- DOC Degradation from either photo- or microbial decomposition
- plant biomass decomposition as % mass loss
TIME POINTS
- Time-0 (T0) = before the addition of plant
materials, plankton were stocked at this stage
- Time-1 (T1) = 10 days after the addition of plant
materials
- Time-2 (T2) = 31 days after …
- Time-3 (T3) = 59 days after …
- Time-4 (T4) = 89 days after …
###DOC
import and for loop to clean up all files and make stacked data in
single df
## import treatment IDs
IDs<-read.csv("data/treatment.IDs.csv")
##### grab files in a list
Total.DOC.files <- list.files(path="data/DOC.TN", pattern = "csv$", full.names = T)
Total.DOC.files## [1] "data/DOC.TN/DOC_T0.csv" "data/DOC.TN/DOC_T1.csv" "data/DOC.TN/DOC_T2.csv"
## [4] "data/DOC.TN/DOC_T3.csv" "data/DOC.TN/DOC_T4.csv"##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="data/DOC.TN", pattern = "csv$", full.names = F))
file.names## [1] "DOC_T0" "DOC_T1" "DOC_T2" "DOC_T3" "DOC_T4"############ formatting all data in for loop
for(i in 1:length(Total.DOC.files))
{
data<-read.csv(Total.DOC.files[i], sep=",")
data<-data[,c(1,3,4)] # removed columns we don't need
data$File<-Total.DOC.files[i]
colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
data$Tank<- IDs$Tank
data$Tank<-as.numeric(as.character(data$Tank)) # make the column of samples all numeric
data <- data[!is.na(as.numeric(as.character(data$Tank))),] # remove all rows that aren't numeric/tanks
data$Treatment<-IDs$Treatment
data$plant.mass..g<-IDs$plant.mass..g
make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
}
########## this is the end of the loop
ls() #see all dfs you've made, the above will be df matching their file names## [1] "data" "DOC_T0" "DOC_T1" "DOC_T2"
## [5] "DOC_T3" "DOC_T4" "Fig.formatting" "file.names"
## [9] "i" "IDs" "Total.DOC.files"#Combine files from loop to single df
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)
DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names
# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27)
#alternatively
# remove the 9 letters ('^.) at start
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File)
# if equals DOC_T0_11052021 then, T0, if not then T1
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
ifelse (DOC.df$File=="DOC_T1.csv", "T1",
ifelse (DOC.df$File=="DOC_T2.csv", "T2",
ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))
#rearrange
DOC.df<- DOC.df %>%
select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L)
DOC.df$Treatment<-as.factor(DOC.df$Treatment)DOC analysis:
- notes on GAMS…. Note that fitting via restricted maximum likelihood (REML) is advised to give stable results. Fitting the smoothing parameter (sp) to user defined levels not advisable (but certainly is done). Number of basis functions important as well, small # of basis functions (k) will limit wiggliness (and make the model more linear than should be). This wiggliness is set by using the k values. Don’t want too high, or too low…
When using the non-linear smoothing, this is the s(xxxx). Can add a new predictor by adding another smooth function with s(x1) + s(x2). When the variable is inside the smooth function, this is to account for the nonlinear shape. Two smoothing together and you have ADDITIVE Nonlinear modeling. Not all variables need to be non-linear, but adding a predictor outside of the s() then this is allowing for linearity.
In practice, continuous variables are rarely linear in GAMs, automatic smoothing will make a linear shape if setting high smoothing function (ie.. s =1000), then a non-linear variable appears linear. However, setting categorical variables as predicots linear is more common.
Factor-smooth interaction also an option, this is when s(x1 by = fac), this sets different smoothing for different variables of “fac”. Usually, the different smoothers are combined with a varying intercept incase the different categories are different in means and slopes, this would be by adding the s(x1 by = fac) + fac, wheere the +fac allows for a different slope
EDF - effective degrees of freedom, equate with wiggliness, with edf =1 as a straight line, and higher edfs as more wiggly. GAM smoother significance is described as not being able to draw a horizontal line through the data. Always check concurvity, which is the collinearity with models from 0-1.
######## T0 model
m1.DOC.T0.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.DOC.T0.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.DOC.T0.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.DOC.AIC<-AIC(m1.DOC.T0.gam, m2.DOC.T0.gam, m3.DOC.T0.gam)
DOC.sum_T0=summary(m3.DOC.T0.gam)
DOC.sum_T0=anova.gam(m3.DOC.T0.gam)
gam.check(m3.DOC.T0.gam, rep=1000)
draw(m3.DOC.T0.gam)
concrvity(m3.DOC.T0.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T0.gam, all.terms = TRUE, page=1)
DOC.diff.T0<-plot_difference(
m1.DOC.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
) +
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T0.mod.plot<-
plot_smooths(
model = m3.DOC.T0.gam,
series = plant.mass..g
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-0") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
######## T1 model
m1.DOC.T1.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.DOC.T1.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.DOC.T1.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.DOC.AIC<-AIC(m1.DOC.T1.gam, m2.DOC.T1.gam, m3.DOC.T1.gam)
summary(m1.DOC.T1.gam)
T1.DOC.anova=anova(m1.DOC.T1.gam)
gam.check(m1.DOC.T1.gam, rep=1000)
draw(m1.DOC.T1.gam)
concrvity(m1.DOC.T1.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T1.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T1<-plot_difference(
m1.DOC.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T1.mod.plot<-
plot_smooths(
model = m1.DOC.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-10") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T2
m1.DOC.T2.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.DOC.T2.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.DOC.T2.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.DOC.AIC<-AIC(m1.DOC.T2.gam, m2.DOC.T2.gam, m3.DOC.T2.gam)
summary(m1.DOC.T2.gam)
T2.DOC.anova=anova.gam(m1.DOC.T2.gam)
gam.check(m1.DOC.T2.gam, rep=1000)
draw(m1.DOC.T2.gam)
concrvity(m1.DOC.T2.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T2.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T2<-plot_difference(
m1.DOC.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T2.mod.plot<-
plot_smooths(
model = m1.DOC.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-31") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T3
m1.DOC.T3.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.DOC.T3.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.DOC.T3.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.DOC.AIC<-AIC(m1.DOC.T3.gam, m2.DOC.T3.gam, m3.DOC.T3.gam)
summary(m1.DOC.T3.gam)
T3.DOC.anova=anova.gam(m1.DOC.T3.gam)
gam.check(m1.DOC.T3.gam, rep=1000)
draw(m1.DOC.T3.gam)
concrvity(m1.DOC.T3.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T3.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T3<-plot_difference(
m1.DOC.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T3.mod.plot<-
plot_smooths(
model = m1.DOC.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-59") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T4
m1.DOC.T4.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.DOC.T4.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.DOC.T4.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.DOC.AIC<-AIC(m1.DOC.T4.gam, m2.DOC.T4.gam, m3.DOC.T4.gam)
summary(m3.DOC.T4.gam)
T4.DOC.anova=anova(m1.DOC.T4.gam)
gam.check(m3.DOC.T4.gam, rep=1000)
draw(m3.DOC.T4.gam)
concrvity(m3.DOC.T4.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T4.gam, all.terms = TRUE, page=1)
###########
#plot for the model output on rawdata
DOC.T4.mod.plot<-
plot_smooths(
model = m3.DOC.T4.gam,
series = plant.mass..g
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-89") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
# model predictions
DOC.diff.T4<-plot_difference(
m1.DOC.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
#labels
mod.DOC.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
#adding columns
mod.DOC.column<- data.frame(mod.DOC.rep)
AIC.DOC<-bind_rows(T1.DOC.AIC, T2.DOC.AIC, T3.DOC.AIC, T4.DOC.AIC)
AIC.DOC.mod<-as.data.frame(cbind(mod.DOC.column, AIC.DOC))
names(AIC.DOC.mod)[1]="Model"
AIC.DOC.mod <- AIC.DOC.mod %>%
# Creating an empty column:
add_column(Variable= c("DOC concentration"), .before="Model")
write.csv(AIC.DOC.mod, "output/AIC models/AIC.DOC.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/DOC.mod.AIC.pdf", width = 6.8, height = 3.8) # Export PDF
grid.table(AIC.DOC.mod, rows = NULL)
dev.off()Generate dataframe for model outputs
####Parametric output
DOC.para.sums=as.data.frame(rbind(T1.DOC.anova$pTerms.table, T2.DOC.anova$pTerms.table, T3.DOC.anova$pTerms.table, T4.DOC.anova$pTerms.table))
names(DOC.para.sums)[1]="df/edf"
#create empty column for ref.df
DOC.para.sums <- DOC.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
DOC.smooth.sums=as.data.frame(rbind(T1.DOC.anova$s.table, T2.DOC.anova$s.table, T3.DOC.anova$s.table, T4.DOC.anova$s.table))
names(DOC.smooth.sums)[1]="df/edf"
#create empty column for ref.df
DOC.smooth.sums <- DOC.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
DOC.mod.allsums=rbind(DOC.para.sums,DOC.smooth.sums)
rownames(DOC.mod.allsums)<-c(1:12)
#create new ID column
DOC.mod.allsums <- DOC.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("DOC concentration"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/DOC.mod.allsums.pdf", width = 9.3, height = 3.8) # Export PDF
grid.table(DOC.mod.allsums, rows= NULL)
dev.off()## quartz_off_screen
## 2DOC: compile raw plots and model-diff plots for final figures
###### compile the plots with effect plots}
extract.legend <- get_legend(
# create some space to the left of the legend
DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))
DOC.alltimes<-plot_grid(
DOC.T0.mod.plot+ theme(legend.position = "none"),
DOC.T1.mod.plot+ theme(legend.position = "none"),
DOC.T2.mod.plot+ theme(legend.position = "none"),
DOC.T3.mod.plot+ theme(legend.position = "none"),
DOC.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/DOC.alltimes.mod.pdf", height=4, width=18)
DOC.alltimes
Figure. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.
DOC.mod.diffs<-plot_grid(
DOC.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
DOC.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
DOC.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
DOC.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
DOC.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8,8), ncol=5)
DOC.mod.diffs
Figure. Model effects from GAMs with differences between smoothers for DOC concentration across time in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.
ggsave("figures/DOC.mod.diffs.pdf", height=4, width=16)UVBIO files: DOC degradation
##### grab files in a list
UVBio.files <- list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = T)
UVBio.files## [1] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T1.csv"
## [2] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T2.csv"
## [3] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T3.csv"
## [4] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T1.csv"
## [5] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T2.csv"
## [6] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T3.csv"##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = F))
file.names## [1] "DOC_NoUV_T1" "DOC_NoUV_T2" "DOC_NoUV_T3" "DOC_UV_T1" "DOC_UV_T2"
## [6] "DOC_UV_T3"############ formatting all data in for loop
for(i in 1:length(UVBio.files))
{
data<-read.csv(UVBio.files[i], sep=",") #skip 13 rows where standards are
data<-data[,c(1,3,4)] # removed columns we don't need
data$File<-UVBio.files[i]
colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
data$Tank<-as.character(data$Tank)
data$Treatment<-IDs$Treatment # add treatment info
data$plant.mass..g<-IDs$plant.mass..g # add plant biomass info
make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
}
########## this is the end of the loop
ls() #see all dfs you've made, the above will be df matching their file names## [1] "AIC.DOC" "AIC.DOC.mod" "data" "DOC_NoUV_T1"
## [5] "DOC_NoUV_T2" "DOC_NoUV_T3" "DOC_T0" "DOC_T1"
## [9] "DOC_T2" "DOC_T3" "DOC_T4" "DOC_UV_T1"
## [13] "DOC_UV_T2" "DOC_UV_T3" "DOC.alltimes" "DOC.df"
## [17] "DOC.diff.T0" "DOC.diff.T1" "DOC.diff.T2" "DOC.diff.T3"
## [21] "DOC.diff.T4" "DOC.mod.allsums" "DOC.mod.diffs" "DOC.para.sums"
## [25] "DOC.smooth.sums" "DOC.sum_T0" "DOC.T0.mod.plot" "DOC.T1.mod.plot"
## [29] "DOC.T2.mod.plot" "DOC.T3.mod.plot" "DOC.T4.mod.plot" "extract.legend"
## [33] "Fig.formatting" "file.names" "i" "IDs"
## [37] "m1.DOC.T0.gam" "m1.DOC.T1.gam" "m1.DOC.T2.gam" "m1.DOC.T3.gam"
## [41] "m1.DOC.T4.gam" "m2.DOC.T0.gam" "m2.DOC.T1.gam" "m2.DOC.T2.gam"
## [45] "m2.DOC.T3.gam" "m2.DOC.T4.gam" "m3.DOC.T0.gam" "m3.DOC.T1.gam"
## [49] "m3.DOC.T2.gam" "m3.DOC.T3.gam" "m3.DOC.T4.gam" "mod.DOC.column"
## [53] "mod.DOC.rep" "T0.DOC.AIC" "T1.DOC.AIC" "T1.DOC.anova"
## [57] "T2.DOC.AIC" "T2.DOC.anova" "T3.DOC.AIC" "T3.DOC.anova"
## [61] "T4.DOC.AIC" "T4.DOC.anova" "Total.DOC.files" "UVBio.files"#Combine files from loop to single df
UVBIO.df<-rbind(DOC_NoUV_T1, DOC_UV_T1, DOC_NoUV_T2, DOC_UV_T2, DOC_NoUV_T3, DOC_UV_T3)
# make new factors for filtered (F/NF) and UV (+/-)
# ** caution! Make sure this matches sheet layout after the loop!
UVBIO.df$Filter.Trt<-rep(c("Filtered", "Unfiltered"), each=30) # repeat 30 times
UVBIO.df$UV.Trt<-rep(c("NoUV", "UV"), each=60) # repeat 60 times
# check tank names in DF and make sure they are what you think they are (i.e 1-30 but w/ weird names)
Tank.fix<-rep(c(1:30), times=12) # repeat 12 times to fix all the strange tank names
#compare tank IDS
df.test<-cbind(Tank.fix, UVBIO.df$Tank) # yep!
UVBIO.df$Tank<-as.factor(Tank.fix) #replace old tank names
# extract out the file names
#Extract file
UVBIO.df$File=sapply(strsplit(UVBIO.df$File, "/"), `[`, 9)
#take out redundancies
UVBIO.df$DOC..mg.L<-gsub("NPOC:","", UVBIO.df$DOC..mg.L)
UVBIO.df$DOC..mg.L<- gsub("mg/L", "", UVBIO.df$DOC..mg.L)
UVBIO.df$TN..mg.L<-gsub("TN:","", UVBIO.df$TN..mg.L)
UVBIO.df$TN..mg.L<- gsub("mg/L", "", UVBIO.df$TN..mg.L)
#convert DOC and TN to numeric
UVBIO.df$DOC..mg.L<-as.numeric(as.character(UVBIO.df$DOC..mg.L))
UVBIO.df$TN..mg.L<- as.numeric(as.character(UVBIO.df$TN..mg.L))
# add a time point
# if equals DOC_UV_T1_111721 then, T1, if not then "nothing"
UVBIO.df$Time.point<- as.factor(ifelse(UVBIO.df$File=="DOC_UV_T1.csv", "T1",
ifelse(UVBIO.df$File=="DOC_NoUV_T1.csv", "T1",
ifelse(UVBIO.df$File=="DOC_UV_T2.csv", "T2",
ifelse(UVBIO.df$File=="DOC_NoUV_T2.csv", "T2", "T3")))))
#rearrange
UVBIO.df<- UVBIO.df %>%
dplyr::select(File, Time.point, Treatment, Filter.Trt, UV.Trt, Tank, plant.mass..g, DOC..mg.L, TN..mg.L)
# make levels for plotting and stats
# add in starting DOC and TN
# original data is... DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4 -- but no need for DOC T0 and T4
# repeat 4x to account for the 4 treatments in each tank
DOC_T1.start<-rep(DOC_T1$DOC..mg.L, times=4)
DOC_T2.start<-rep(DOC_T2$DOC..mg.L, times=4)
DOC_T3.start<-rep(DOC_T3$DOC..mg.L, times=4)
TN_T1.start<-rep(DOC_T1$TN..mg.L, times=4)
TN_T2.start<-rep(DOC_T2$TN..mg.L, times=4)
TN_T3.start<-rep(DOC_T3$TN..mg.L, times=4)
# compile and add in as a new column
UVBIO.df$Start.DOC<-c(DOC_T1.start, DOC_T2.start, DOC_T3.start)
UVBIO.df$Start.TN<-c(TN_T1.start, TN_T2.start, TN_T3.start)
UVBIO.df$Burn.Filt.UV<-interaction(UVBIO.df$Treatment, UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)
UVBIO.df$Filt.UV<-interaction(UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)###SUVA SUVA calculations Time point 0
#load file and rearrange
suva.t0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T0.csv")
suva.t0.df$Sample.Name=gsub("\\..*", "", suva.t0.df$Sample.Name)
names(suva.t0.df)[1]="Tank"
suva.t0.df=suva.t0.df[-1,]
suva.t0.df$Tank=gsub("T", "", suva.t0.df$Tank)
suva.t0.df$Tank=as.numeric(suva.t0.df$Tank)
suva.t0.df=suva.t0.df[order(suva.t0.df$Tank),]
#add important info
suva.t0.df$DOC..mg.L=DOC_T0$DOC..mg.L
suva.t0.df$Time.point="T0"
suva.t0.df$Treatment=IDs$Treatment
suva.t0.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t0.df$SUVA=(suva.t0.df[suva.t0.df$Tank,]$abs254/suva.t0.df[suva.t0.df$Tank,]$DOC..mg.L)*100Time point 1
#load file and rearrange
library(dplyr)
suva.t1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T1.csv")
suva.t1.df=suva.t1.df[-1,]
suva.t1.df=suva.t1.df[-1,]
names(suva.t1.df)[1]="Tank"
suva.t1.df$Tank=gsub("TA", "", suva.t1.df$Tank)
suva.t1.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t1.df$Tank)
suva.t1.df$Tank=sub("\\..*", "", suva.t1.df$Tank)
suva.t1.df$Tank=as.numeric(suva.t1.df$Tank)
suva.t1.df=suva.t1.df[order(suva.t1.df$Tank),]
#add important info
suva.t1.df$DOC..mg.L=DOC_T1$DOC..mg.L
suva.t1.df$Time.point="T1"
suva.t1.df$Treatment=IDs$Treatment
suva.t1.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t1.df$SUVA=(suva.t1.df[suva.t1.df$Tank,]$abs254/suva.t1.df[suva.t1.df$Tank,]$DOC..mg.L)*100Time point 2
#load file and rearrange
suva.t2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T2_d1A.edit.csv")
names(suva.t2.df)[1]="Tank"
suva.t2.df$Tank=gsub("TA", "", suva.t2.df$Tank)
suva.t2.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t2.df$Tank)
suva.t2.df$Tank=sub("\\..*", "", suva.t2.df$Tank)
suva.t2.df$Tank=as.numeric(suva.t2.df$Tank)
suva.t2.df=suva.t2.df[order(suva.t2.df$Tank),]
#add important info
suva.t2.df$DOC..mg.L=DOC_T2$DOC..mg.L
suva.t2.df$Time.point="T2"
suva.t2.df$Treatment=IDs$Treatment
suva.t2.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t2.df$SUVA=(suva.t2.df[suva.t2.df$Tank,]$abs254/suva.t2.df[suva.t2.df$Tank,]$DOC..mg.L)*100Time point 3
suva.t3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T3.csv")
names(suva.t3.df)[1]="Tank"
suva.t3.df$Tank=gsub("TA", "", suva.t3.df$Tank)
suva.t3.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t3.df$Tank)
suva.t3.df$Tank=sub("\\..*", "", suva.t3.df$Tank)
suva.t3.df$Tank=as.numeric(suva.t3.df$Tank)
suva.t3.df=suva.t3.df[order(suva.t3.df$Tank),]
#add important info
suva.t3.df$DOC..mg.L=DOC_T3$DOC..mg.L
suva.t3.df$Time.point="T3"
suva.t3.df$Treatment=IDs$Treatment
suva.t3.df$plant.mass..g=IDs$plant.mass..g
class(suva.t3.df$DOC..mg.L)## [1] "numeric"#calculate SUVA
suva.t3.df$SUVA=(suva.t3.df[suva.t3.df$Tank,]$abs254/suva.t3.df[suva.t3.df$Tank,]$DOC..mg.L)*100Time point 4
suva.t4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T4.csv")
names(suva.t4.df)[1]="Tank"
suva.t4.df$Tank=gsub("TA", "", suva.t4.df$Tank)
suva.t4.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t4.df$Tank)
suva.t4.df$Tank=sub("\\..*", "", suva.t4.df$Tank)
suva.t4.df$Tank=as.numeric(suva.t4.df$Tank)
suva.t4.df=suva.t4.df[order(suva.t4.df$Tank),]
#add important info
suva.t4.df$DOC..mg.L=DOC_T4$DOC..mg.L
suva.t4.df$Time.point="T4"
suva.t4.df$Treatment=IDs$Treatment
suva.t4.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t4.df$SUVA=(suva.t4.df[suva.t4.df$Tank,]$abs254/suva.t4.df[suva.t4.df$Tank,]$DOC..mg.L)*100Combine all SUVA time points
suva.df=rbind(suva.t0.df, suva.t1.df, suva.t2.df, suva.t3.df, suva.t4.df)Other indices
# Time 0
index.T0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T0_indices.csv")
index.T0.df$Time.point="T0"
names(index.T0.df)[6]="Humification.Index"
index.T0.df=index.T0.df[-c(7:11,13,14)]
#Time 1
index.T1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T1_indices.csv")
index.T1.df$Time.point="T1"
index.T1.df=index.T1.df[-c(7:11,13,14)]
#Time 2
index.T2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T2_indices.csv")
index.T2.df$Time.point="T2"
index.T2.df=index.T2.df[-c(7:11,13,14)]
#Time 3
index.T3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T3_indices.csv")
index.T3.df$Time.point="T3"
index.T3.df=index.T3.df[-c(7:11,13,14)]
#Time 4
index.T4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T4_indices.csv")
index.T4.df$Time.point="T4"
index.T4.df=index.T4.df[-c(7:11,13,14)]
#combine all index time points
index.df=rbind(index.T0.df, index.T1.df, index.T2.df, index.T3.df, index.T4.df)Add these new metrics to DOC df import and do a loop to clean up all files and make stacked data in single df
#Combine files from loop to single df
DOC_T0$SUVA=suva.df[suva.df$Time.point=="T0",]$SUVA
DOC_T1$SUVA=suva.df[suva.df$Time.point=="T1",]$SUVA
DOC_T2$SUVA=suva.df[suva.df$Time.point=="T2",]$SUVA
DOC_T3$SUVA=suva.df[suva.df$Time.point=="T3",]$SUVA
DOC_T4$SUVA=suva.df[suva.df$Time.point=="T4",]$SUVA
#HIX
DOC_T0$hix=index.df[index.df$Time.point=="T0",]$Humification.Index
DOC_T1$hix=index.df[index.df$Time.point=="T1",]$Humification.Index
DOC_T2$hix=index.df[index.df$Time.point=="T2",]$Humification.Index
DOC_T3$hix=index.df[index.df$Time.point=="T3",]$Humification.Index
DOC_T4$hix=index.df[index.df$Time.point=="T4",]$Humification.Index
#Fluorescence
DOC_T0$fl=index.df[index.df$Time.point=="T0",]$Fluroscence.Index
DOC_T1$fl=index.df[index.df$Time.point=="T1",]$Fluroscence.Index
DOC_T2$fl=index.df[index.df$Time.point=="T2",]$Fluroscence.Index
DOC_T3$fl=index.df[index.df$Time.point=="T3",]$Fluroscence.Index
DOC_T4$fl=index.df[index.df$Time.point=="T4",]$Fluroscence.Index
#slope ratio
DOC_T0$sr=suva.df[suva.df$Time.point=="T0",]$Spectral.Slope.Ratio
DOC_T1$sr=suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio
DOC_T2$sr=suva.df[suva.df$Time.point=="T2",]$Spectral.Slope.Ratio
DOC_T3$sr=suva.df[suva.df$Time.point=="T3",]$Spectral.Slope.Ratio
DOC_T4$sr=suva.df[suva.df$Time.point=="T4",]$Spectral.Slope.Ratio
#freshness (BIX)
DOC_T0$fr=index.df[index.df$Time.point=="T0",]$Freshness.Index
DOC_T1$fr=index.df[index.df$Time.point=="T1",]$Freshness.Index
DOC_T2$fr=index.df[index.df$Time.point=="T2",]$Freshness.Index
DOC_T3$fr=index.df[index.df$Time.point=="T3",]$Freshness.Index
DOC_T4$fr=index.df[index.df$Time.point=="T4",]$Freshness.Index
#combine
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)
#Extract file
DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names
# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27)
#alternatively
# remove the 9 letters ('^.) at start
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File)
# if equals DOC_T0 then, T0, if not then T1..etc
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
ifelse (DOC.df$File=="DOC_T1.csv", "T1",
ifelse (DOC.df$File=="DOC_T2.csv", "T2",
ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))
#rearrange
library(dplyr)
DOC.df<- DOC.df %>%
dplyr::select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L, SUVA, sr, hix, fl, fr)
DOC.df=DOC.df%>%
mutate_if(is.character,as.factor)SUVA analysis and plots
######## T0 model
m1.SUVA.T0.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.SUVA.T0.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.SUVA.T0.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.suva.AIC<-AIC(m1.SUVA.T0.gam, m2.SUVA.T0.gam, m3.SUVA.T0.gam)
T0.suva.AIC## df AIC
## m1.SUVA.T0.gam 8.544088 68.24603
## m2.SUVA.T0.gam 4.730649 70.94187
## m3.SUVA.T0.gam 3.000364 81.69088summary(m1.SUVA.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3007 0.1664 13.824 1.04e-12 ***
## Treatmentunburned 1.0159 0.2354 4.316 0.000251 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.000 0.995 0.329
## s(plant.mass..g):Treatmentunburned 3.734 4.544 1.598 0.176
##
## R-sq.(adj) = 0.46 Deviance explained = 56.7%
## -REML = 34.112 Scale est. = 0.41545 n = 30anova.gam(m1.SUVA.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric Terms:
## df F p-value
## Treatment 1 18.63 0.000251
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.000 0.995 0.329
## s(plant.mass..g):Treatmentunburned 3.734 4.544 1.598 0.176gam.check(m1.SUVA.T0.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.233511e-05,6.231324e-06]
## (score 34.11213 & scale 0.4154513).
## Hessian positive definite, eigenvalue range [1.233492e-05,13.1468].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.00 1.06 0.49
## s(plant.mass..g):Treatmentunburned 9.00 3.73 1.06 0.56draw(m1.SUVA.T0.gam)
concrvity(m1.SUVA.T0.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 1.91e-29
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.SUVA.T0.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T0<-plot_difference(
m1.SUVA.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
suva.T0.mod.plot<-
plot_smooths(
model = m3.SUVA.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
coord_cartesian(ylim = c(0,12))+
xlab("plant material (g)") +
Fig.formatting
suva.T0.mod.plot
######## T1 model
m1.suva.T1.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.suva.T1.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.suva.T1.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.suva.AIC<-AIC(m1.suva.T1.gam, m2.suva.T1.gam, m3.suva.T1.gam)
T1.suva.AIC## df AIC
## m1.suva.T1.gam 11.829790 16.68484
## m2.suva.T1.gam 7.758699 18.60925
## m3.suva.T1.gam 6.762722 17.68119summary(m1.suva.T1.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.28324 0.06804 33.559 <2e-16 ***
## Treatmentunburned 0.09731 0.09622 1.011 0.324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 4.597 5.541 20.70 < 2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.354 4.095 14.88 6.98e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.857 Deviance explained = 90.1%
## -REML = 15.005 Scale est. = 0.069436 n = 30T1.suva.anova=anova(m1.suva.T1.gam)
gam.check(m1.suva.T1.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-6.613102e-06,3.156923e-06]
## (score 15.00512 & scale 0.06943626).
## Hessian positive definite, eigenvalue range [0.9040944,13.38198].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 4.60 0.77 0.025 *
## s(plant.mass..g):Treatmentunburned 9.00 3.35 0.77 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1draw(m1.suva.T1.gam)
concrvity(m1.suva.T1.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 1.18e-29
## 6 observed s(plant.mass..g):Treatmentunburned 2.90e-30
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.suva.T1.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T1<-plot_difference(
m1.suva.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T1.mod.plot<-
plot_smooths(
model = m1.suva.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
coord_cartesian(ylim = c(0,12))+
xlab("") +
Fig.formatting
########## T2
m1.suva.T2.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.suva.T2.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.suva.T2.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.suva.AIC<-AIC(m1.suva.T2.gam, m2.suva.T2.gam, m3.suva.T2.gam)
T2.suva.AIC## df AIC
## m1.suva.T2.gam 6.560914 166.8072
## m2.suva.T2.gam 5.095612 168.0063
## m3.suva.T2.gam 4.087052 167.3995summary(m1.suva.T2.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9241 0.8879 5.546 9.24e-06 ***
## Treatmentunburned -1.4665 1.2557 -1.168 0.254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.075 2.561 2.763 0.065 .
## s(plant.mass..g):Treatmentunburned 1.000 1.000 0.068 0.796
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.17 Deviance explained = 28.7%
## -REML = 75.261 Scale est. = 11.826 n = 30T2.suva.anova=anova(m1.suva.T2.gam)
gam.check(m1.suva.T2.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-3.815864e-05,5.391568e-06]
## (score 75.26059 & scale 11.82582).
## Hessian positive definite, eigenvalue range [3.815639e-05,13.02274].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.07 1.11 0.65
## s(plant.mass..g):Treatmentunburned 9.00 1.00 1.11 0.56draw(m1.suva.T2.gam)
concrvity(m1.suva.T2.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 1.54e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.84e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.suva.T2.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T2<-plot_difference(
m1.suva.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T2.mod.plot<-
plot_smooths(
model = m1.suva.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0, 12) +
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
suva.T2.mod.plot
########## T3
m1.suva.T3.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.suva.T3.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.suva.T3.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.suva.AIC<-AIC(m1.suva.T3.gam, m2.suva.T3.gam, m3.suva.T3.gam)
T3.suva.AIC## df AIC
## m1.suva.T3.gam 10.018360 23.87907
## m2.suva.T3.gam 7.067145 28.31202
## m3.suva.T3.gam 5.518278 36.92467summary(m1.suva.T3.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.95603 0.07731 38.238 < 2e-16 ***
## Treatmentunburned 0.39851 0.10933 3.645 0.00141 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.458 3.023 76.33 <2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.270 3.996 92.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.955 Deviance explained = 96.5%
## -REML = 14.721 Scale est. = 0.089645 n = 30T3.suva.anova=anova(m1.suva.T3.gam)
gam.check(m1.suva.T3.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [3.454876e-10,1.6981e-08]
## (score 14.72117 & scale 0.0896449).
## Hessian positive definite, eigenvalue range [0.2269231,13.14457].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.46 0.81 0.14
## s(plant.mass..g):Treatmentunburned 9.00 3.27 0.81 0.11draw(m1.suva.T3.gam)
concrvity(m1.suva.T3.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 4.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 1.39e-30
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.suva.T3.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T3<-plot_difference(
m1.suva.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T3.mod.plot<-
plot_smooths(
model = m1.suva.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,12))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.suva.T4.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.suva.T4.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.suva.T4.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.suva.AIC<-AIC(m1.suva.T4.gam, m2.suva.T4.gam, m3.suva.T4.gam)
T4.suva.AIC## df AIC
## m1.suva.T4.gam 6.949281 33.36689
## m2.suva.T4.gam 5.560335 38.15182
## m3.suva.T4.gam 4.576257 36.69035# best is global
summary(m1.suva.T4.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.66735 0.09569 27.874 <2e-16 ***
## Treatmentunburned -0.10136 0.13533 -0.749 0.461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.00 130.62 <2e-16 ***
## s(plant.mass..g):Treatmentunburned 2.531 3.11 52.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.909 Deviance explained = 92.3%
## -REML = 17.897 Scale est. = 0.13736 n = 30T4.suva.anova=anova(m1.suva.T4.gam)
gam.check(m1.suva.T4.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-8.722995e-06,1.789474e-06]
## (score 17.89694 & scale 0.1373614).
## Hessian positive definite, eigenvalue range [8.722846e-06,13.04785].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.00 0.98 0.35
## s(plant.mass..g):Treatmentunburned 9.00 2.53 0.98 0.34draw(m1.suva.T4.gam)
concrvity(m1.suva.T4.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 8.97e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.suva.T4.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T4<-plot_difference(
m1.suva.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T4.mod.plot<-
plot_smooths(
model = m1.suva.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,12))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
suva.T4.mod.plot
mod.suva.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.suva.column<- data.frame(mod.suva.rep)
AIC.suva<-bind_rows(T1.suva.AIC, T2.suva.AIC, T3.suva.AIC, T4.suva.AIC)
AIC.suva.mod<-as.data.frame(cbind(mod.suva.column, AIC.suva))
names(AIC.suva.mod)[1]="Model"
AIC.suva.mod <- AIC.suva.mod %>%
# Creating an empty column:
add_column(Variable= c("SUVA"), .before="Model")
write.csv(AIC.suva.mod, "output/AIC models/AIC.suva.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/suva.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.suva.mod, rows= NULL)
dev.off()## quartz_off_screen
## 2Generate dataframe for model outputs
####Parametric output
suva.para.sums=as.data.frame(rbind(T1.suva.anova$pTerms.table, T2.suva.anova$pTerms.table, T3.suva.anova$pTerms.table, T4.suva.anova$pTerms.table))
names(suva.para.sums)[1]="df/edf"
suva.para.sums
#create empty column for ref.df
suva.para.sums <- suva.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
suva.smooth.sums=as.data.frame(rbind(T1.suva.anova$s.table, T2.suva.anova$s.table, T3.suva.anova$s.table, T4.suva.anova$s.table))
names(suva.smooth.sums)[1]="df/edf"
#create empty column for ref.df
suva.smooth.sums <- suva.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
suva.mod.allsums=rbind(suva.para.sums,suva.smooth.sums)
rownames(suva.mod.allsums)<-c(1:12)
#create new ID column
suva.mod.allsums <- suva.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("SUVA"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/suva.mod.allsums.pdf", width = 8.5, height = 3.8) # Export PDF
grid.table(suva.mod.allsums, rows= NULL)
dev.off()compile raw plots and model-diff plots for final figures
extract.legend <- get_legend(
# create some space to the left of the legend
DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))
suva.alltimes<-plot_grid(
suva.T1.mod.plot+ theme(legend.position = "none"),
suva.T2.mod.plot+ theme(legend.position = "none"),
suva.T3.mod.plot+ theme(legend.position = "none"),
suva.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
suva.alltimes
ggsave("figures/SUVA.alltimes.mod.pdf", height=4, width=12)
suva.mod.diffs<-plot_grid(
suva.diff.T1+ theme(legend.position = "none")+ ggtitle("SUVA - Day-10"),
suva.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
suva.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
suva.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
suva.mod.diffs
ggsave("figures/suva.mod.diffs.pdf", height=4, width=12)Sr analysis and plots
######## T0 model
m1.sr.T0.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.sr.T0.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.sr.T0.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.sr.AIC<-AIC(m1.sr.T0.gam, m2.sr.T0.gam, m3.sr.T0.gam)
T0.sr.AIC## df AIC
## m1.sr.T0.gam 5.656040 72.39828
## m2.sr.T0.gam 4.000309 76.34838
## m3.sr.T0.gam 3.000189 77.53584summary(m1.sr.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8530 0.1871 20.591 <2e-16 ***
## Treatmentunburned -0.5062 0.2646 -1.913 0.067 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.000 6.896 0.0143 *
## s(plant.mass..g):Treatmentunburned 1.376 1.656 0.401 0.5296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.234 Deviance explained = 32.3%
## -REML = 34.155 Scale est. = 0.52519 n = 30anova.gam(m1.sr.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric Terms:
## df F p-value
## Treatment 1 3.659 0.067
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.000 6.896 0.0143
## s(plant.mass..g):Treatmentunburned 1.376 1.656 0.401 0.5296gam.check(m1.sr.T0.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.375453e-05,2.556971e-06]
## (score 34.15521 & scale 0.5251928).
## Hessian positive definite, eigenvalue range [1.375414e-05,13.00273].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.00 0.9 0.19
## s(plant.mass..g):Treatmentunburned 9.00 1.38 0.9 0.24draw(m1.sr.T0.gam)
concrvity(m1.sr.T0.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.79e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.sr.T0.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T0<-plot_difference(
m1.sr.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
sr.T0.mod.plot<-
plot_smooths(
model = m3.sr.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylab(expression('S'[italic(R)]))+
xlab("plant material (g)") +
Fig.formatting
sr.T0.mod.plot
######## T1 model
m1.sr.T1.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.sr.T1.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.sr.T1.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.sr.AIC<-AIC(m1.sr.T1.gam, m2.sr.T1.gam, m3.sr.T1.gam)
T1.sr.AIC## df AIC
## m1.sr.T1.gam 7.127378 7.775483
## m2.sr.T1.gam 5.172277 12.562474
## m3.sr.T1.gam 4.186417 11.180366summary(m1.sr.T1.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.29977 0.06211 20.93 <2e-16 ***
## Treatmentunburned 0.07200 0.08784 0.82 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.000 1.000 12.85 0.00149 **
## s(plant.mass..g):Treatmentunburned 2.545 3.127 13.36 2.25e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.643 Deviance explained = 69.9%
## -REML = 6.6791 Scale est. = 0.057874 n = 30T1.sr.anova=anova(m1.sr.T1.gam)
gam.check(m1.sr.T1.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-2.139521e-06,2.098115e-07]
## (score 6.679058 & scale 0.05787405).
## Hessian positive definite, eigenvalue range [2.139513e-06,13.04761].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.00 0.65 0.015 *
## s(plant.mass..g):Treatmentunburned 9.00 2.55 0.65 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1draw(m1.sr.T1.gam)
concrvity(m1.sr.T1.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 7.24e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.sr.T1.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T1<-plot_difference(
m1.sr.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T1.mod.plot<-
plot_smooths(
model = m1.sr.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab(expression('S'[italic(R)]))+
coord_cartesian(ylim = c(0.25,2.5))+
xlab("") +
Fig.formatting
sr.T1.mod.plot
########## T2
m1.sr.T2.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.sr.T2.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.sr.T2.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.sr.AIC<-AIC(m1.sr.T2.gam, m2.sr.T2.gam, m3.sr.T2.gam)
T2.sr.AIC## df AIC
## m1.sr.T2.gam 11.537024 -13.51026
## m2.sr.T2.gam 7.677441 -20.11997
## m3.sr.T2.gam 6.712723 -22.00780summary(m3.sr.T2.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ s(plant.mass..g)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4040 0.0269 52.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g) 4.171 5.053 17.25 8.53e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.748 Deviance explained = 78.4%
## -REML = -6.5014 Scale est. = 0.021714 n = 30T2.sr.anova=anova(m1.sr.T2.gam)
gam.check(m3.sr.T2.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.027443e-09,5.61081e-10]
## (score -6.501377 & scale 0.02171413).
## Hessian positive definite, eigenvalue range [1.070181,14.19444].
## Model rank = 10 / 10
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g) 9.00 4.17 1.59 1draw(m3.sr.T2.gam)
concrvity(m3.sr.T2.gam)## # A tibble: 6 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 4.61e-25
## 2 worst s(plant.mass..g) 4.61e-25
## 3 observed para 4.61e-25
## 4 observed s(plant.mass..g) 3.94e-31
## 5 estimate para 4.61e-25
## 6 estimate s(plant.mass..g) 1.24e-27par(mfrow = c(2, 2))
plot(m3.sr.T2.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T2<-plot_difference(
m1.sr.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T2.mod.plot<-
plot_smooths(
model = m3.sr.T2.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.25,2.5))+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
sr.T2.mod.plot 
########## T3
m1.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.sr.T3.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m4.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "ML")
T3.sr.AIC<-AIC(m1.sr.T3.gam, m2.sr.T3.gam, m3.sr.T3.gam)
T3.sr.AIC## df AIC
## m1.sr.T3.gam 9.688689 -1.124062
## m2.sr.T3.gam 6.776426 -6.818488
## m3.sr.T3.gam 5.747507 -6.468986summary(m1.sr.T3.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.39164 0.05136 27.098 <2e-16 ***
## Treatmentunburned -0.09869 0.07263 -1.359 0.188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 3.035 3.715 11.99 4.06e-05 ***
## s(plant.mass..g):Treatmentunburned 2.548 3.131 11.13 0.000108 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.729 Deviance explained = 79%
## -REML = 3.8444 Scale est. = 0.039562 n = 30T3.sr.anova=anova(m1.sr.T3.gam)
gam.check(m1.sr.T3.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-1.513811e-11,1.06315e-12]
## (score 3.844409 & scale 0.03956209).
## Hessian positive definite, eigenvalue range [0.6287135,13.13229].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 3.04 1.05 0.52
## s(plant.mass..g):Treatmentunburned 9.00 2.55 1.05 0.50draw(m1.sr.T3.gam)
concrvity(m1.sr.T3.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 1.64e-30
## 6 observed s(plant.mass..g):Treatmentunburned 6.91e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.sr.T3.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T3<-plot_difference(
m1.sr.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T3.mod.plot<-
plot_smooths(
model = m2.sr.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen"))+
coord_cartesian(ylim = c(0.25,2.5))+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
sr.T3.mod.plot
########## T4
m1.sr.T4.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.sr.T4.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.sr.T4.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.sr.AIC<-AIC(m1.sr.T4.gam, m2.sr.T4.gam, m3.sr.T4.gam)
T4.sr.AIC## df AIC
## m1.sr.T4.gam 11.398371 -57.07811
## m2.sr.T4.gam 7.363534 -42.13080
## m3.sr.T4.gam 6.396737 -44.12181summary(m1.sr.T4.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.288325 0.019910 64.709 <2e-16 ***
## Treatmentunburned 0.004562 0.028156 0.162 0.873
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.890 3.542 25.64 <2e-16 ***
## s(plant.mass..g):Treatmentunburned 4.498 5.427 46.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.922 Deviance explained = 94.4%
## -REML = -17.839 Scale est. = 0.0059459 n = 30T4.sr.anova=anova(m1.sr.T4.gam)
gam.check(m1.sr.T4.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-5.197224e-09,4.086935e-10]
## (score -17.83887 & scale 0.005945903).
## Hessian positive definite, eigenvalue range [0.6347134,13.32902].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.89 1 0.45
## s(plant.mass..g):Treatmentunburned 9.00 4.50 1 0.47draw(m1.sr.T4.gam)
concrvity(m1.sr.T4.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 3.17e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.19e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.sr.T4.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T4<-plot_difference(
m1.sr.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T4.mod.plot<-
plot_smooths(
model = m1.sr.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.25,2.5))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
sr.T4.mod.plot
mod.sr.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.sr.column<- data.frame(mod.sr.rep)
AIC.sr<-bind_rows(T1.sr.AIC, T2.sr.AIC, T3.sr.AIC, T4.sr.AIC)
AIC.sr## df AIC
## m1.sr.T1.gam 7.127378 7.775483
## m2.sr.T1.gam 5.172277 12.562474
## m3.sr.T1.gam 4.186417 11.180366
## m1.sr.T2.gam 11.537024 -13.510260
## m2.sr.T2.gam 7.677441 -20.119965
## m3.sr.T2.gam 6.712723 -22.007797
## m1.sr.T3.gam 9.688689 -1.124062
## m2.sr.T3.gam 6.776426 -6.818488
## m3.sr.T3.gam 5.747507 -6.468986
## m1.sr.T4.gam 11.398371 -57.078114
## m2.sr.T4.gam 7.363534 -42.130799
## m3.sr.T4.gam 6.396737 -44.121805AIC.sr.mod<-as.data.frame(cbind(mod.sr.column, AIC.sr))
names(AIC.sr.mod)[1]="Model"
AIC.sr.mod <- AIC.sr.mod %>%
# Creating an empty column:
add_column(Variable= c("Slope ratio"), .before="Model")
write.csv(AIC.sr.mod, "output/AIC models/AIC.sr.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/sr.mod.AIC.pdf", width = 6.3, height = 3.8) # Export PDF
grid.table(AIC.sr.mod, rows=NULL)
dev.off()## quartz_off_screen
## 2Generate dataframe for model outputs
####Parametric output
sr.para.sums=as.data.frame(rbind(T1.sr.anova$pTerms.table, T2.sr.anova$pTerms.table, T3.sr.anova$pTerms.table, T4.sr.anova$pTerms.table))
names(sr.para.sums)[1]="df/edf"
#create empty column for ref.df
sr.para.sums <- sr.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
sr.smooth.sums=as.data.frame(rbind(T1.sr.anova$s.table, T2.sr.anova$s.table, T3.sr.anova$s.table, T4.sr.anova$s.table))
names(sr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
sr.smooth.sums <- sr.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
sr.smooth.sums
#merge the dataframes
sr.mod.allsums=rbind(sr.para.sums,sr.smooth.sums)
rownames(sr.mod.allsums)<-c(1:12)
#create new ID column
sr.mod.allsums <- sr.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("Slope ratio"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/sr.mod.allsums.pdf", width = 8.5, height = 3.8) # Export PDF
grid.table(sr.mod.allsums, rows= NULL)
dev.off()Sr: compile raw plots and model-diff plots for final figures
sr.alltimes<-plot_grid(
sr.T1.mod.plot+ theme(legend.position = "none"),
sr.T2.mod.plot+ theme(legend.position = "none"),
sr.T3.mod.plot+ theme(legend.position = "none"),
sr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
sr.alltimes
ggsave("figures/sr.alltimes.mod.pdf", height=4, width=12)
sr.mod.diffs<-plot_grid(
sr.diff.T1+ theme(legend.position = "none")+ ggtitle("Sr - Day-10"),
sr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
sr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
sr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
sr.mod.diffs
ggsave("figures/sr.mod.diffs.pdf", height=4, width=12)HIX analysis and plots
######## T0 model
m1.hix.T0.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.hix.T0.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.hix.T0.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.hix.AIC<-AIC(m1.hix.T0.gam, m2.hix.T0.gam, m3.hix.T0.gam)
T0.hix.AIC## df AIC
## m1.hix.T0.gam 5.700361 20.37557
## m2.hix.T0.gam 5.240725 18.48445
## m3.hix.T0.gam 3.985728 23.07642summary(m1.hix.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.52625 0.07855 19.429 <2e-16 ***
## Treatmentunburned -0.27629 0.11109 -2.487 0.0197 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.406 1.7 0.489 0.693
## s(plant.mass..g):Treatmentunburned 1.000 1.0 2.258 0.145
##
## R-sq.(adj) = 0.178 Deviance explained = 27.4%
## -REML = 11.609 Scale est. = 0.092562 n = 30anova.gam(m1.hix.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric Terms:
## df F p-value
## Treatment 1 6.185 0.0197
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.406 1.700 0.489 0.693
## s(plant.mass..g):Treatmentunburned 1.000 1.000 2.258 0.145gam.check(m1.hix.T0.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-6.467201e-06,1.475016e-05]
## (score 11.60873 & scale 0.09256152).
## Hessian positive definite, eigenvalue range [6.46737e-06,13.00317].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.41 0.95 0.39
## s(plant.mass..g):Treatmentunburned 9.00 1.00 0.95 0.28draw(m1.hix.T0.gam)
concrvity(m1.hix.T0.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 9.13e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.84e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.hix.T0.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T0<-plot_difference(
m1.hix.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
hix.T0.mod.plot<-
plot_smooths(
model = m3.hix.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
coord_cartesian(ylim = c(0,6))+
ylab("HIX")+
xlab("plant material (g)") +
Fig.formatting
######## T1 model
m1.hix.T1.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.hix.T1.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.hix.T1.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.hix.AIC<-AIC(m1.hix.T1.gam, m2.hix.T1.gam, m3.hix.T1.gam)
T1.hix.AIC## df AIC
## m1.hix.T1.gam 5.000061 62.10694
## m2.hix.T1.gam 4.000060 60.23250
## m3.hix.T1.gam 3.000061 59.42114summary(m3.hix.T1.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ s(plant.mass..g)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7544 0.1114 15.75 1.92e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g) 1 1 15.64 0.000475 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.335 Deviance explained = 35.8%
## -REML = 29.297 Scale est. = 0.37225 n = 30T1.hix.anova=anova(m1.hix.T1.gam)
gam.check(m3.sr.T1.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-6.202238e-07,-3.062326e-07]
## (score 6.592856 & scale 0.07082567).
## Hessian positive definite, eigenvalue range [0.2360873,14.01068].
## Model rank = 10 / 10
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g) 9.00 1.77 0.73 0.035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1draw(m3.hix.T1.gam)
concrvity(m3.hix.T1.gam)## # A tibble: 6 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 4.61e-25
## 2 worst s(plant.mass..g) 4.61e-25
## 3 observed para 4.61e-25
## 4 observed s(plant.mass..g) 2.60e-32
## 5 estimate para 4.61e-25
## 6 estimate s(plant.mass..g) 1.24e-27par(mfrow = c(2, 2))
plot(m3.hix.T1.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T1<-plot_difference(
m1.hix.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T1.mod.plot<-
plot_smooths(
model = m3.hix.T1.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-10") +
ylab("HIX")+
xlab("") +
Fig.formatting
hix.T1.mod.plot
########## T2
m1.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.hix.T2.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m4.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.hix.AIC<-AIC(m1.hix.T2.gam, m2.hix.T2.gam, m3.hix.T2.gam)
T2.hix.AIC## df AIC
## m1.hix.T2.gam 5.000022 30.35861
## m2.hix.T2.gam 4.000030 28.41308
## m3.hix.T2.gam 3.000045 28.88814summary(m2.hix.T2.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ Treatment + s(plant.mass..g)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.17247 0.09254 23.475 <2e-16 ***
## Treatmentunburned 0.19943 0.13088 1.524 0.139
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g) 1 1 6.651 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.194 Deviance explained = 24.9%
## -REML = 15.017 Scale est. = 0.12847 n = 30T2.hix.anova=anova(m1.hix.T2.gam)
gam.check(m2.hix.T2.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 13 iterations.
## Gradient range [-5.861103e-06,1.733725e-06]
## (score 15.01688 & scale 0.1284671).
## Hessian positive definite, eigenvalue range [5.861019e-06,13.5].
## Model rank = 11 / 11
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g) 9 1 0.86 0.12draw(m2.hix.T2.gam)
concrvity(m2.hix.T2.gam)## # A tibble: 6 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g) 4.61e-25
## 3 observed para 5 e- 1
## 4 observed s(plant.mass..g) 2.96e-32
## 5 estimate para 5 e- 1
## 6 estimate s(plant.mass..g) 1.24e-27par(mfrow = c(2, 2))
plot(m2.hix.T2.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T2<-plot_difference(
m1.hix.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T2.mod.plot<-
plot_smooths(
model = m2.hix.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
hix.T2.mod.plot
########## T3
m1.hix.T3.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.hix.T3.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.hix.T3.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.hix.AIC<-AIC(m1.hix.T3.gam, m2.hix.T3.gam, m3.hix.T3.gam)
T3.hix.AIC## df AIC
## m1.hix.T3.gam 11.953543 23.11911
## m2.hix.T3.gam 7.757426 23.50220
## m3.hix.T3.gam 5.936863 35.52011summary(m1.hix.T3.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4393 0.0745 32.741 < 2e-16 ***
## Treatmentunburned 0.4249 0.1054 4.033 0.000623 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.841 3.484 7.129 0.001285 **
## s(plant.mass..g):Treatmentunburned 4.609 5.553 6.346 0.000791 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.718 Deviance explained = 80%
## -REML = 16.611 Scale est. = 0.083257 n = 30T3.hix.anova=anova(m1.hix.T3.gam)
gam.check(m1.hix.T3.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-1.41462e-08,1.860842e-09]
## (score 16.61062 & scale 0.08325724).
## Hessian positive definite, eigenvalue range [0.6660674,13.3327].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.84 0.91 0.26
## s(plant.mass..g):Treatmentunburned 9.00 4.61 0.91 0.23draw(m1.hix.T3.gam)
concrvity(m1.hix.T3.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 9.86e-31
## 6 observed s(plant.mass..g):Treatmentunburned 2.60e-30
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.hix.T3.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T3<-plot_difference(
m1.hix.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T3.mod.plot<-
plot_smooths(
model = m1.hix.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.hix.T4.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.hix.T4.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.hix.T4.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.hix.AIC<-AIC(m1.hix.T4.gam, m2.hix.T4.gam, m3.hix.T4.gam)
T4.hix.AIC## df AIC
## m1.hix.T4.gam 10.053318 25.12183
## m2.hix.T4.gam 6.208475 41.21670
## m3.hix.T4.gam 4.841730 52.43251summary(m1.hix.T4.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.59052 0.07888 32.840 < 2e-16 ***
## Treatmentunburned 0.58969 0.11156 5.286 2.55e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.501 3.074 10.75 0.000141 ***
## s(plant.mass..g):Treatmentunburned 3.257 3.980 37.18 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.877 Deviance explained = 90.5%
## -REML = 15.28 Scale est. = 0.09334 n = 30T4.hix.anova=anova(m1.hix.T4.gam)
gam.check(m1.hix.T4.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-6.675986e-06,2.083557e-06]
## (score 15.28009 & scale 0.0933395).
## Hessian positive definite, eigenvalue range [0.3536202,13.14571].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.50 0.73 0.070 .
## s(plant.mass..g):Treatmentunburned 9.00 3.26 0.73 0.045 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1draw(m1.hix.T4.gam)
concrvity(m1.hix.T4.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 1.95e-30
## 6 observed s(plant.mass..g):Treatmentunburned 2.31e-30
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.hix.T4.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T4<-plot_difference(
m1.hix.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T4.mod.plot<-
plot_smooths(
model = m1.hix.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
mod.hix.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.hix.column<- data.frame(mod.hix.rep)
AIC.hix<-bind_rows(T1.hix.AIC, T2.hix.AIC, T3.hix.AIC, T4.hix.AIC)
AIC.hix.mod<-as.data.frame(cbind(mod.hix.column, AIC.hix))
names(AIC.hix.mod)[1]="Model"
AIC.hix.mod <- AIC.hix.mod %>%
# Creating an empty column:
add_column(Variable= c("HIX"), .before="Model")
write.csv(AIC.hix.mod, "output/AIC models/AIC.hix.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/hix.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.hix.mod, rows= NULL)
dev.off()## quartz_off_screen
## 2Generate dataframe for model outputs
####Parametric output
hix.para.sums=as.data.frame(rbind(T1.hix.anova$pTerms.table, T2.hix.anova$pTerms.table, T3.hix.anova$pTerms.table, T4.hix.anova$pTerms.table))
names(hix.para.sums)[1]="df/edf"
#create empty column for ref.df
hix.para.sums <- hix.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
hix.smooth.sums=as.data.frame(rbind(T1.hix.anova$s.table, T2.hix.anova$s.table, T3.hix.anova$s.table, T4.hix.anova$s.table))
names(hix.smooth.sums)[1]="df/edf"
#create empty column for ref.df
hix.smooth.sums <- hix.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
hix.mod.allsums=rbind(hix.para.sums,hix.smooth.sums)
#create new ID column
hix.mod.allsums <- hix.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("HIX"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/hix.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(hix.mod.allsums, rows= NULL)
dev.off()## quartz_off_screen
## 2HIX: compile raw plots and model-diff plots for final figures
hix.alltimes<-plot_grid(
hix.T1.mod.plot+ theme(legend.position = "none"),
hix.T2.mod.plot+ theme(legend.position = "none"),
hix.T3.mod.plot+ theme(legend.position = "none"),
hix.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
hix.alltimes
ggsave("figures/hix.alltimes.mod.pdf", height=4, width=12)
hix.mod.diffs<-plot_grid(
hix.diff.T1+ theme(legend.position = "none")+ ggtitle("HIX - Day-10"),
hix.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
hix.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
hix.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
hix.mod.diffs
ggsave("figures/hix.mod.diffs.pdf", height=4, width=12)Freshness analysis and plots
######## T0 model
m1.fr.T0.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.fr.T0.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.fr.T0.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.fr.AIC<-AIC(m1.fr.T0.gam, m2.fr.T0.gam, m3.fr.T0.gam)
T0.fr.AIC## df AIC
## m1.fr.T0.gam 5.000114 -56.31473
## m2.fr.T0.gam 4.000035 -56.48287
## m3.fr.T0.gam 3.000036 -58.25380summary(m1.fr.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.89630 0.02222 40.33 <2e-16 ***
## Treatmentunburned 0.01446 0.03143 0.46 0.649
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1 1 1.038 0.318
## s(plant.mass..g):Treatmentunburned 1 1 0.625 0.436
##
## R-sq.(adj) = -0.0404 Deviance explained = 6.73%
## -REML = -21.46 Scale est. = 0.0074074 n = 30anova.gam(m1.fr.T0.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric Terms:
## df F p-value
## Treatment 1 0.212 0.649
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1 1 1.038 0.318
## s(plant.mass..g):Treatmentunburned 1 1 0.625 0.436gam.check(m1.fr.T0.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-1.342278e-05,1.476944e-05]
## (score -21.46003 & scale 0.007407422).
## Hessian positive definite, eigenvalue range [4.232612e-07,12.99999].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9 1 1.13 0.75
## s(plant.mass..g):Treatmentunburned 9 1 1.13 0.72draw(m1.fr.T0.gam)
concrvity(m1.fr.T0.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned 5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.84e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.fr.T0.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T0<-plot_difference(
m1.fr.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
fr.T0.mod.plot<-
plot_smooths(
model = m3.fr.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylim(0.25,1.75)+
xlab("plant material (g)") +
Fig.formatting
fr.T0.mod.plot
######## T1 model
m1.fr.T1.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.fr.T1.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.fr.T1.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.fr.AIC<-AIC(m1.fr.T1.gam, m2.fr.T1.gam, m3.fr.T1.gam)
T1.fr.AIC## df AIC
## m1.fr.T1.gam 6.194913 -5.4939710
## m2.fr.T1.gam 4.000010 -0.6819142
## m3.fr.T1.gam 3.000009 2.7456950summary(m1.fr.T1.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.67640 0.05057 13.375 5.89e-13 ***
## Treatmentunburned 0.18648 0.07152 2.607 0.0151 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 1.774 2.195 0.971 0.446767
## s(plant.mass..g):Treatmentunburned 1.000 1.000 18.832 0.000207 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.46 Deviance explained = 53%
## -REML = 0.45251 Scale est. = 0.038362 n = 30T1.fr.anova=anova(m1.fr.T1.gam)
gam.check(m1.fr.T1.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-8.295506e-07,2.641631e-06]
## (score 0.4525131 & scale 0.0383616).
## Hessian positive definite, eigenvalue range [8.295565e-07,13.01166].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 1.77 0.9 0.22
## s(plant.mass..g):Treatmentunburned 9.00 1.00 0.9 0.14draw(m1.fr.T1.gam)
concrvity(m1.fr.T1.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 6.67e-31
## 6 observed s(plant.mass..g):Treatmentunburned 9.84e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.fr.T1.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T1<-plot_difference(
m1.fr.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T1.mod.plot<-
plot_smooths(
model = m1.fr.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylim(0.25,1.75)+
ylab(expression("Freshness"~"("~italic(alpha)~":"~italic(beta)~")"))+
xlab("") +
Fig.formatting
fr.T1.mod.plot
########## T2
m1.fr.T2.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.fr.T2.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.fr.T2.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.fr.AIC<-AIC(m1.fr.T2.gam, m2.fr.T2.gam, m3.fr.T2.gam)
T2.fr.AIC## df AIC
## m1.fr.T2.gam 6.756227 -102.78590
## m2.fr.T2.gam 5.460197 -101.45652
## m3.fr.T2.gam 4.296831 -96.00129summary(m1.fr.T2.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.790989 0.009899 79.907 < 2e-16 ***
## Treatmentunburned -0.039854 0.013999 -2.847 0.00874 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g):Treatmentburned 2.237 2.756 3.361 0.03249 *
## s(plant.mass..g):Treatmentunburned 1.000 1.000 13.326 0.00121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.493 Deviance explained = 56.7%
## -REML = -41.46 Scale est. = 0.0014698 n = 30T2.fr.anova=anova(m1.fr.T2.gam)
gam.check(m1.fr.T2.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.292444e-05,6.966117e-06]
## (score -41.46041 & scale 0.001469797).
## Hessian positive definite, eigenvalue range [1.292412e-05,13.03067].
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g):Treatmentburned 9.00 2.24 1.06 0.48
## s(plant.mass..g):Treatmentunburned 9.00 1.00 1.06 0.57draw(m1.fr.T2.gam)
concrvity(m1.fr.T2.gam)## # A tibble: 9 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g):Treatmentburned 4.60e-25
## 3 worst s(plant.mass..g):Treatmentunburned 4.58e-25
## 4 observed para 5 e- 1
## 5 observed s(plant.mass..g):Treatmentburned 8.83e-30
## 6 observed s(plant.mass..g):Treatmentunburned 9.84e-31
## 7 estimate para 5 e- 1
## 8 estimate s(plant.mass..g):Treatmentburned 1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned 1.24e-27par(mfrow = c(2, 2))
plot(m1.fr.T2.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T2<-plot_difference(
m1.fr.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T2.mod.plot<-
plot_smooths(
model = m1.fr.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
fr.T2.mod.plot
########## T3
m1.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.fr.T3.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m4.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.fr.AIC<-AIC(m1.fr.T3.gam, m2.fr.T3.gam, m3.fr.T3.gam)
T3.fr.AIC## df AIC
## m1.fr.T3.gam 8.337961 -104.51469
## m2.fr.T3.gam 6.186728 -108.46478
## m3.fr.T3.gam 4.617005 -86.80427summary(m2.fr.T3.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.815763 0.009067 89.969 < 2e-16 ***
## Treatmentunburned -0.071997 0.012823 -5.615 7.29e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g) 2.631 3.231 8.329 0.000457 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.661 Deviance explained = 70.3%
## -REML = -46.175 Scale est. = 0.0012332 n = 30T3.fr.anova=anova(m1.fr.T3.gam)
gam.check(m2.fr.T3.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.323279e-09,8.61311e-11]
## (score -46.1752 & scale 0.001233207).
## Hessian positive definite, eigenvalue range [0.6071085,13.5516].
## Model rank = 11 / 11
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g) 9.00 2.63 1.26 0.91draw(m2.fr.T3.gam)
concrvity(m2.fr.T3.gam)## # A tibble: 6 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g) 4.61e-25
## 3 observed para 5 e- 1
## 4 observed s(plant.mass..g) 2.12e-31
## 5 estimate para 5 e- 1
## 6 estimate s(plant.mass..g) 1.24e-27par(mfrow = c(2, 2))
plot(m2.fr.T3.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T3<-plot_difference(
m1.fr.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T3.mod.plot<-
plot_smooths(
model = m2.fr.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.fr.T4.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m4.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.fr.AIC<-AIC(m1.fr.T4.gam, m2.fr.T4.gam, m3.fr.T4.gam)
T4.fr.AIC## df AIC
## m1.fr.T4.gam 9.278312 -103.54444
## m2.fr.T4.gam 7.825392 -104.36395
## m3.fr.T4.gam 3.470877 -93.80125summary(m2.fr.T4.gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## fr ~ Treatment + s(plant.mass..g)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.814924 0.009446 86.27 < 2e-16 ***
## Treatmentunburned -0.043554 0.013359 -3.26 0.00331 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant.mass..g) 3.975 4.825 7.147 0.000483 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.593 Deviance explained = 66.3%
## -REML = -42.963 Scale est. = 0.0013384 n = 30T4.fr.anova=anova(m1.fr.T4.gam)
gam.check(m2.fr.T4.gam, rep=1000)
##
## Method: REML Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-3.584282e-06,1.105561e-06]
## (score -42.96326 & scale 0.001338422).
## Hessian positive definite, eigenvalue range [0.2725441,13.66728].
## Model rank = 11 / 11
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(plant.mass..g) 9.00 3.98 1.08 0.56draw(m2.fr.T4.gam)
concrvity(m2.fr.T4.gam)## # A tibble: 6 × 3
## type term concurvity
## <chr> <chr> <dbl>
## 1 worst para 5 e- 1
## 2 worst s(plant.mass..g) 4.61e-25
## 3 observed para 5 e- 1
## 4 observed s(plant.mass..g) 4.41e-31
## 5 estimate para 5 e- 1
## 6 estimate s(plant.mass..g) 1.24e-27par(mfrow = c(2, 2))
plot(m2.fr.T4.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T4<-plot_difference(
m1.fr.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T4.mod.plot<-
plot_smooths(
model = m2.fr.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
ggtitle("Day-89") +
geom_line(aes(linetype=Treatment))+
ylab("") +
xlab("") +
Fig.formattingGenerate dataframe for model outputs
####Parametric output
fr.para.sums=as.data.frame(rbind(T1.fr.anova$pTerms.table, T2.fr.anova$pTerms.table, T3.fr.anova$pTerms.table, T4.fr.anova$pTerms.table))
names(fr.para.sums)[1]="df/edf"
#create empty column for ref.df
fr.para.sums <- fr.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
fr.smooth.sums=as.data.frame(rbind(T1.fr.anova$s.table, T2.fr.anova$s.table, T3.fr.anova$s.table, T4.fr.anova$s.table))
names(fr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fr.smooth.sums <- fr.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
fr.mod.allsums=rbind(fr.para.sums,fr.smooth.sums)
#create new ID column
fr.mod.allsums <- fr.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("FR"), .before="Model")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(fr.mod.allsums, rows= NULL)
dev.off()
mod.fr.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.fr.column<- data.frame(mod.fr.rep)
AIC.fr<-bind_rows(T1.fr.AIC, T2.fr.AIC, T3.fr.AIC, T4.fr.AIC)
AIC.fr.mod<-as.data.frame(cbind(mod.fr.column, AIC.fr))
names(AIC.fr.mod)[1]="Model"
AIC.fr.mod <- AIC.fr.mod %>%
# Creating an empty column:
add_column(Variable= c("FR"), .before="Model")
write.csv(AIC.fr.mod, "output/AIC models/AIC.fr.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.AIC.pdf", width = 6.3, height = 3.8) # Export PDF
grid.table(AIC.fr.mod, rows= NULL)
dev.off()Freshness: compile raw plots and model-diff plots for final figures
fr.alltimes<-plot_grid(
fr.T1.mod.plot+ theme(legend.position = "none"),
fr.T2.mod.plot+ theme(legend.position = "none"),
fr.T3.mod.plot+ theme(legend.position = "none"),
fr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
fr.alltimes
ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)
fr.mod.diffs<-plot_grid(
fr.diff.T1+ theme(legend.position = "none")+ ggtitle("FR - Day-10"),
fr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
fr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
fr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
fr.mod.diffs
ggsave("figures/fr.mod.diffs.pdf", height=4, width=12)Fluorescence analysis and plots
library(mgcv)
library(tidymv)
######## T0 model
m1.fl.T0.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.fl.T0.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.fl.T0.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.fl.AIC<-AIC(m1.fl.T0.gam, m2.fl.T0.gam, m3.fl.T0.gam)
T0.fl.AIC
summary(m1.fl.T0.gam)
anova.gam(m1.fl.T0.gam)
gam.check(m1.fl.T0.gam, rep=1000)
draw(m1.fl.T0.gam)
concrvity(m1.fl.T0.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T0.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T0<-plot_difference(
m1.fl.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,0.5))+
Fig.formatting
###########
fl.T0.mod.plot<-
plot_smooths(
model = m1.fl.T0.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
coord_cartesian(ylim = c(0.5,2.5))+
ylab("FI")+
xlab("plant material (g)") +
Fig.formatting
fl.T0.mod.plot
######## T1 model
m1.fl.T1.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.fl.T1.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.fl.T1.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.fl.AIC<-AIC(m1.fl.T1.gam, m2.fl.T1.gam, m3.fl.T1.gam)
T1.fl.AIC
summary(m1.fl.T1.gam)
T1.fl.anova=anova(m1.fl.T1.gam)
T1.fl.anova
gam.check(m1.fl.T1.gam, rep=1000)
draw(m1.fl.T1.gam)
concrvity(m1.fl.T1.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T1.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T1<-plot_difference(
m1.fl.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T1.mod.plot<-
plot_smooths(
model = m1.fl.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab("FI")+
xlab("") +
ylim(0.5,2.5)+
Fig.formatting
fl.T1.mod.plot
########## T2
m1.fl.T2.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.fl.T2.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.fl.T2.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.fl.AIC<-AIC(m1.fl.T2.gam, m2.fl.T2.gam, m3.fl.T2.gam)
T2.fl.AIC
summary(m1.fl.T2.gam)
T2.fl.anova=anova(m1.fl.T2.gam)
gam.check(m1.fl.T2.gam, rep=1000)
draw(m1.fl.T2.gam)
concrvity(m1.fl.T2.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T2.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T2<-plot_difference(
m1.fl.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T2.mod.plot<-
plot_smooths(
model = m1.fl.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
coord_cartesian(ylim = c(0.5,2.5))+
Fig.formatting
fl.T2.mod.plot
########## T3
m1.fl.T3.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.fl.T3.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.fl.T3.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.fl.AIC<-AIC(m1.fl.T3.gam, m2.fl.T3.gam, m3.fl.T3.gam)
T3.fl.AIC
summary(m3.fl.T3.gam)
T3.fl.anova=anova(m1.fl.T3.gam)
gam.check(m3.fl.T3.gam, rep=1000)
draw(m3.fl.T3.gam)
concrvity(m3.fl.T3.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T3.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T3<-plot_difference(
m1.fl.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T3.mod.plot<-
plot_smooths(
model = m3.fl.T3.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.fl.T4.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.fl.T4.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.fl.T4.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.fl.AIC<-AIC(m1.fl.T4.gam, m2.fl.T4.gam, m3.fl.T4.gam)
T4.fl.AIC
summary(m3.fl.T4.gam)
T4.fl.anova=anova(m1.fl.T4.gam)
gam.check(m3.fl.T4.gam, rep=1000)
draw(m3.fl.T4.gam)
concrvity(m3.fl.T4.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T4.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T4<-plot_difference(
m1.fl.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T4.mod.plot<-
plot_smooths(
model = m3.fl.T4.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
mod.fl.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.fl.column<- data.frame(mod.fl.rep)
AIC.fl<-bind_rows(T1.fl.AIC, T2.fl.AIC, T3.fl.AIC, T4.fl.AIC)
AIC.fl.mod<-as.data.frame(cbind(mod.fl.column, AIC.fl))
AIC.fl.mod
names(AIC.fl.mod)[1]="Model"
AIC.fl.mod <- AIC.fl.mod %>%
# Creating an empty column:
add_column(Variable= c("FI"), .before="Model")
write.csv(AIC.fl.mod, "output/AIC models/AIC.fl.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.fl.mod, rows=NULL)
dev.off()Generate dataframe for model outputs
####Parametric output
fl.para.sums=as.data.frame(rbind(T1.fl.anova$pTerms.table, T2.fl.anova$pTerms.table, T3.fl.anova$pTerms.table, T4.fl.anova$pTerms.table))
names(fl.para.sums)[1]="df/edf"
#create empty column for ref.df
fl.para.sums <- fl.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
fl.smooth.sums=as.data.frame(rbind(T1.fl.anova$s.table, T2.fl.anova$s.table, T3.fl.anova$s.table, T4.fl.anova$s.table))
names(fl.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fl.smooth.sums <- fl.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
fl.mod.allsums=rbind(fl.para.sums,fl.smooth.sums)
rownames(fl.mod.allsums)<-c(1:12)
#create new ID column
fl.mod.allsums <- fl.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("FI"), .before="Model")
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(fl.mod.allsums, rows=NULL)
dev.off()Fluorescence: compile raw plots and model-diff plots for final figures
fl.alltimes<-plot_grid(
fl.T1.mod.plot+ theme(legend.position = "none"),
fl.T2.mod.plot+ theme(legend.position = "none"),
fl.T3.mod.plot+ theme(legend.position = "none"),
fl.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
fl.alltimes
ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)
fl.mod.diffs<-plot_grid(
fl.diff.T1+ theme(legend.position = "none")+ ggtitle("FI - Day-10"),
fl.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
fl.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
fl.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
fl.mod.diffs
ggsave("figures/fl.mod.diffs.pdf", height=4, width=12)Compile all indices into one plot
indices.alltimes= plot_grid(suva.alltimes, sr.alltimes, hix.alltimes, fr.alltimes, fl.alltimes, nrow=5, rel_widths = c(8,8,8,8,8), labels = "AUTO")
ggsave("figures/indices.alltimes.pdf", height=18, width=14)
indices.diff.plot= plot_grid(suva.mod.diffs, sr.mod.diffs, hix.mod.diffs, fr.mod.diffs, fl.mod.diffs, nrow=5, rel_widths = c(8,8,8,8,8), labels = "auto")
indices.diff.plot
ggsave("figures/indices.diff.plot.pdf", height=18, width=14)
#### Day 0 baseline indicies
day0_index.baseline.plot = plot_grid(
suva.T0.mod.plot+ theme(legend.position = "none"),
sr.T0.mod.plot+ theme(legend.position = "none"),
hix.T0.mod.plot+ theme(legend.position = "none"),
fr.T0.mod.plot+ theme(legend.position = "none"),
fl.T0.mod.plot+ theme(legend.position = "none"),
extract.legend,
rel_widths = c(8,8,8,8,8,3), nrow=6)
day0_index.baseline.plot
ggsave("figures/day0_index.baseline.plot.pdf", height=18, width=14)
Calculate photo, microbe, and total degradation:
#+UV/+Bio=A
#-UV/+Bio=B
#+UV/-Bio=C
#-UV/-Bio (NoUV,Filtered)=D
#This is calculating % decrease ((control-treatment)/control)*100
#photodegradation
photo.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L),
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
photo.df$effect="UV-only"
names(photo.df)[1]="Time.point"
names(photo.df)[2]="Start.DOC"
names(photo.df)[3]="Start.TN"
names(photo.df)[4]="Treatment"
names(photo.df)[5]="Tank"
names(photo.df)[6]="plant.mass..g"
names(photo.df)[7]="abs.change.DOC"
names(photo.df)[8]="percent.change.DOC"
names(photo.df)[9]="abs.change.TN"
names(photo.df)[10]="percent.change.TN"
names(photo.df)[11]="SUVA"
names(photo.df)[12]="Spectral.slope"
names(photo.df)[13]="HIX"
names(photo.df)[14]="Fr"
names(photo.df)[15]="FI"
photo.df$mass.remaining.DOC=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.DOC)
photo.df$mass.remaining.TN=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.TN)
photo_t1.df=subset(photo.df, Time.point %in% c("T1"))
photo.df=photo.df%>%
mutate_if(is.character,as.factor)
photo_t1.df=photo_t1.df%>%
mutate_if(is.character,as.factor)
#Microbe decomp
microbe.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L- UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(( UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
microbe.df$effect="microbes-only"
names(microbe.df)[1]="Time.point"
names(microbe.df)[2]="Start.DOC"
names(microbe.df)[3]="Start.TN"
names(microbe.df)[4]="Treatment"
names(microbe.df)[5]="Tank"
names(microbe.df)[6]="plant.mass..g"
names(microbe.df)[7]="abs.change.DOC"
names(microbe.df)[8]="percent.change.DOC"
names(microbe.df)[9]="abs.change.TN"
names(microbe.df)[10]="percent.change.TN"
names(microbe.df)[11]="SUVA"
names(microbe.df)[12]="Spectral.slope"
names(microbe.df)[13]="HIX"
names(microbe.df)[14]="Fr"
names(microbe.df)[15]="FI"
microbe.df$mass.remaining.DOC=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.DOC)
microbe.df$mass.remaining.TN=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.TN)
microbe_t1.df=subset(microbe.df, Time.point %in% c("T1"))
microbe.df=microbe.df%>%
mutate_if(is.character,as.factor)
microbe_t1.df=microbe_t1.df%>%
mutate_if(is.character,as.factor)
#Total decomp
total.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
total.df$effect="UV and microbes"
names(total.df)[1]="Time.point"
names(total.df)[2]="Start.DOC"
names(total.df)[3]="Start.TN"
names(total.df)[4]="Treatment"
names(total.df)[5]="Tank"
names(total.df)[6]="plant.mass..g"
names(total.df)[7]="abs.change.DOC"
names(total.df)[8]="percent.change.DOC"
names(total.df)[9]="abs.change.TN"
names(total.df)[10]="percent.change.TN"
names(total.df)[11]="SUVA"
names(total.df)[12]="Spectral.slope"
names(total.df)[13]="HIX"
names(total.df)[14]="Fr"
names(total.df)[15]="FI"
total.df$mass.remaining.DOC=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.DOC)
total.df$mass.remaining.TN=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.TN)
total_t1.df=subset(total.df, Time.point %in% c("T1"))
total.df=total.df%>%
mutate_if(is.character,as.factor)
total_t1.df=total_t1.df%>%
mutate_if(is.character,as.factor)
#merge degra.df time points into one dataframe
degra.df=rbind(photo.df, microbe.df, total.df)
degra.df=pivot_longer(degra.df, cols=c("effect"), names_to="effect", values_to="Type")
degra_t1.df=rbind(photo_t1.df, microbe_t1.df, total_t1.df)
degra_t1.df=pivot_longer(degra_t1.df, cols=c("effect"), names_to="effect", values_to="Type")
#make sure we have factors
degra.df=degra.df%>%
mutate_if(is.character,as.factor)
degra_t1.df=degra_t1.df%>%
mutate_if(is.character,as.factor)
#seperate degra.df for burned and unburned tanks
degra_t1.burned.df=degra_t1.df[degra_t1.df$Treatment=="burned",]
degra_t1.unburned.df=degra_t1.df[degra_t1.df$Treatment=="unburned",]UVBIO models and plots
library(mgcv)
###Just UV
m1.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m2.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m3.photo.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m1.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m2.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m3.photo.abs=gam(abs.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
photo.AIC=AIC(m1.photo.perc, m2.photo.perc, m3.photo.perc)
photo.AIC=AIC(m1.photo.abs, m2.photo.abs, m3.photo.abs)
photo.AIC
summary(m1.photo.perc)
photo.anova=anova(m1.photo.perc)
photo.anova
gam.check(m3.photo.perc, rep=1000)
draw(m3.photo.perc)
concrvity(m3.photo.perc)
par(mfrow = c(2, 2))
plot(m3.photo.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
photo.mod.plot<-
plot_smooths(
model = m1.photo.perc,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=photo.df[(photo.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
geom_hline(yintercept=0, linetype="longdash", color = "gray")+
ggtitle("UV-only") +
ylab("Change DOC (%)") +
xlab("plant material (g)") +
coord_cartesian(ylim = c(-30,10))+
Fig.formatting
photo.mod.plot
###Just microbes
m1.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
m2.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
m3.microbe.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
microbe.AIC=AIC(m1.microbe.perc, m2.microbe.perc, m3.microbe.perc)
microbe.AIC
summary(m3.microbe.perc)
microbe.anova=anova(m1.microbe.perc)
gam.check(m3.microbe.perc, rep=1000)
draw(m3.microbe.perc)
concrvity(m3.microbe.perc)
par(mfrow = c(2, 2))
plot(m3.microbe.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
microbe.mod.plot<-
plot_smooths(
model = m3.microbe.perc,
series = plant.mass..g
) +
geom_point(data=microbe.df[(microbe.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("microbes-only") +
geom_hline(yintercept=0, linetype="longdash", color = "gray") +
ylab("Change DOC (%)") +
coord_cartesian(ylim = c(-30,10))+
xlab("plant material (g)") +
Fig.formatting
coord_cartesian(ylim = c(-30,0))
###Just total decomp
m1.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m2.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m3.total.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m4.total.perc=gam(percent.change.DOC ~ s(plant.mass..g, m=c(2,0)), data=total.df, method="REML", select= T, subset= Time.point=="T1")
total.AIC=AIC(m1.total.perc, m2.total.perc, m3.total.perc)
total.AIC
summary(m3.total.perc)
total.anova=anova(m1.total.perc)
gam.check(m3.total.perc, rep=1000)
draw(m3.total.perc)
concrvity(m3.total.perc)
par(mfrow = c(2, 2))
plot(m3.total.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
total.mod.plot<-
plot_smooths(
model = m3.total.perc,
series = plant.mass..g
) +
geom_point(data=total.df[(total.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("UV and microbes") +
geom_hline(yintercept=0, linetype="longdash", color = "gray") +
ylab("Change DOC (%)") +
coord_cartesian(ylim = c(-30,10))+
xlab("plant material (g)") +
Fig.formatting
##combine AIC
mod.degra.rep<-rep(c("UV-only: Treatment with by-factor smooth",
"UV-only: Treatment with global smooth",
"UV-only: global smooth",
"Microbes-only: Treatment with by-factor smooth",
"Microbes-only: Treatment with global smooth",
"Microbes-only: global smooth",
"UV and microbes: Treatment with by-factor smooth",
"UV and microbes: Treatment with global smooth",
"UV and microbes: global smooth"))
mod.degra.column<- data.frame(mod.degra.rep)
AIC.degra<-bind_rows(photo.AIC, microbe.AIC, total.AIC)
AIC.degra.mod<-as.data.frame(cbind(mod.degra.column, AIC.degra))
AIC.degra.mod
names(AIC.degra.mod)[1]="Model"
AIC.degra.mod <- AIC.degra.mod %>%
# Creating an empty column:
add_column(Variable= c("-"), .before="Model")
write.csv(AIC.degra.mod, "output/AIC models/AIC.degra.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.AIC.pdf", width = 6.5, height = 2.8) # Export PDF
grid.table(AIC.degra.mod, rows=NULL)
dev.off()
calculate effect size for degradation pathways
library(effectsize)
r_squared_microbes <- summary(m1.microbe.perc)$r.sq
r_squared_total <- summary(m1.total.perc)$r.sq
r_squared_microbes## [1] 0.3219554r_squared_total## [1] 0.4629789Generate dataframe for model outputs
####Parametric output
degra.para.sums=as.data.frame(rbind(photo.anova$pTerms.table, microbe.anova$pTerms.table, total.anova$pTerms.table))
names(degra.para.sums)[1]="df/edf"
#create empty column for ref.df
degra.para.sums <- degra.para.sums %>%
# Creating an empty column:
add_column(Model= c("UV-only Model: Treatment","Microbes-only Model: Treatment","UV and microbes Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
degra.smooth.sums=as.data.frame(rbind(photo.anova$s.table, microbe.anova$s.table, total.anova$s.table))
names(degra.smooth.sums)[1]="df/edf"
#create empty column for ref.df
degra.smooth.sums <- degra.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("UV-only Model: burned smooth","UV-only Model: unburned smooth","Microbes-only Model: burned smooth","Microbes-only Model: unburned smooth","UV and microbes Model: burned smooth","UV and microbes Model: unburned smooth"), .before="df/edf")
#merge the dataframes
degra.mod.allsums=rbind(degra.para.sums,degra.smooth.sums)
degra.mod.allsums <- degra.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("-"), .before="Model")
write.csv(degra.mod.allsums,"/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.pdf", width = 8.5, height = 2.8) # Export PDF
grid.table(degra.mod.allsums, rows=NULL)
dev.off()degra mod plots
degra.mod.plots = plot_grid(
photo.mod.plot + theme(legend.position = "none"),
microbe.mod.plot + theme(legend.position = "none"),
total.mod.plot + theme(legend.position = "none" ),
extract.legend,
rel_widths = c(8,8,8,3), nrow = 1)
degra.mod.plots
ggsave("figures/degra.mod.plots.pdf", height=4, width=10)Sage and Willow Decomposition Analysis
#load file
plant.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania2/data/decomp.weights.csv")
plant.df=plant.df[-c(61:82), ]
plant.df=plant.df[-c(6, 8,9)]
#make sure we have factors
plant.df=plant.df%>%
mutate_if(is.character,as.factor)
#add propotion column so it can be arcsine transformed
plant.df$perc.loss=(100*(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g-plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$dry.mass..g.bag.corrected)/(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g))
#remove NaN
plant.df[is.na(plant.df)] <- 0
#calculate mass remaining
plant.df$mass.remaining=100-plant.df$perc.loss
hist(plant.df$mass.remaining)
#add delta
plant.df$delta.loss= plant.df$perc.loss/100
#arcsine
plant.df$delta.trans=asin(sqrt(plant.df$delta.loss))
hist(plant.df$delta.trans)
#exclude control tank outlier
plant.df2=data.frame(plant.df[plant.df$plant.specific.mass..g > 0,])
###Sage
m1.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m2.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m3.sage.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m4.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
hist(m3.sage.gam$residuals)
m5.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m6.sage.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
sage.AIC=AIC(m4.sage.gam, m5.sage.gam, m6.sage.gam)
sage.AIC
sage.anova=anova(m4.sage.gam)
sage.anova
gam.check(m6.sage.gam, rep=1000)
mean(plant.df[plant.df$plant.type=="sage",]$delta.trans)
# global smooth
# model predictions
sage.diff.mod<-plot_difference(
m1.sage.gam,
series = plant.specific.mass..g,
difference = list(treatment = c("burned", "unburned"))
)
sage.diff.mod
summary(m4.sage.gam)
anova.gam(m3.sage.gam)
gam.check(m4.sage.gam, rep=1000)
draw(m3.sage.gam)
concrvity(m3.sage.gam)
par(mfrow = c(2, 2))
plot(m3.sage.gam, all.terms = TRUE, page=1)
qq_plot(m3.sage.gam)
#plot for the model output on rawdata
sage.mod.plot<-
plot_smooths(
model = m4.sage.gam,
series = plant.specific.mass..g,
comparison= treatment
) +
geom_point(data=plant.df2[(plant.df2$plant.type=="sage"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Sage") +
ylab("Mass loss (%)") +
coord_cartesian(ylim = c(0,80))+
xlab("Plant material (g)") +
Fig.formatting
sage.mod.plot
### Willow
m1.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m2.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m3.willow.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m4.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m5.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df, method="REML", select= T, subset= plant.type=="willow")
m6.willow.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
willow.AIC=AIC( m4.willow.gam,m5.willow.gam,m6.willow.gam)
willow.AIC
# treatment with global smooth
summary(m4.willow.gam)
willow.anova=anova(m4.willow.gam)
willow.anova
gam.check(m3.willow.gam, rep=1000)
draw(m2.willow.gam)
concrvity(m2.willow.gam)
par(mfrow = c(2, 2))
plot(m2.willow.gam, all.terms = TRUE, page=1)
# model predictions
willow.diff.mod<-plot_difference(
m1.willow.gam,
series = plant.specific.mass..g,
difference = list(treatment = c("burned", "unburned"))
)
willow.diff.mod
#plot for the model output on rawdata
willow.mod.plot<-
plot_smooths(
model = m4.willow.gam,
series = plant.specific.mass..g,
comparison = treatment
) +
geom_point(data=plant.df2[(plant.df2$plant.type=="willow"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Willow") +
ylab("Mass loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
willow.mod.plot
##combine AIC
mod.decomp.rep<-rep(c("Sage: Treatment with by-factor smooth",
"Sage: Treatment with global smooth",
"Sage: global smooth",
"Willow: Treatment with by-factor smooth",
"Willow: Treatment with global smooth",
"Willow: global smooth"))
mod.decomp.column<- data.frame(mod.decomp.rep)
AIC.decomp<-bind_rows(sage.AIC, willow.AIC)
AIC.decomp.mod<-as.data.frame(cbind(mod.decomp.column, AIC.decomp))
names(AIC.decomp.mod)[1]="Model"
write.csv(AIC.decomp.mod, "output/AIC models/AIC.decomp.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.AIC.pdf", width = 5, height = 2) # Export PDF
grid.table(AIC.decomp.mod, rows=NULL)
dev.off()
Generate dataframe for model outputs
####Parametric output
decomp.para.sums=as.data.frame(rbind(sage.anova$pTerms.table, willow.anova$pTerms.table))
names(decomp.para.sums)[1]="df/edf"
#create empty column for ref.df
decomp.para.sums <- decomp.para.sums %>%
# Creating an empty column:
add_column(Model= c("Sage Model: Treatment","Willow Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
decomp.smooth.sums=as.data.frame(rbind(sage.anova$s.table, willow.anova$s.table))
names(decomp.smooth.sums)[1]="df/edf"
#create empty column for ref.df
decomp.smooth.sums <- decomp.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Sage Model: burned smooth","Sage Model: unburned smooth","Willow Model: burned smooth","Willow Model: unburned smooth"), .before="df/edf")
#merge the dataframes
decomp.mod.allsums=rbind(decomp.para.sums,decomp.smooth.sums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.allsums.pdf", width = 7, height = 2) # Export PDF
grid.table(decomp.mod.allsums, rows=NULL)
dev.off()
comparing sage and willow by fire treatment
#compare mass loss between sage and willow by fire treatment---burned
m1.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
m2.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
m3.compare.burned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
compare.burned.AIC=AIC(m1.compare.burned.gam, m2.compare.burned.gam, m3.compare.burned.gam)
compare.burned.anova=anova(m1.compare.burned.gam)
compare.burned.anova
gam.check(m1.compare.burned.gam, rep=1000)
# model predictions
compare.burned.diff.mod<-plot_difference(
m1.compare.burned.gam,
series = plant.specific.mass..g,
difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod
#plot for the model output on rawdata
compare.burned.mod.plot<-
plot_smooths(
model = m2.compare.burned.gam,
series = plant.specific.mass..g,
comparison = plant.type
) +
geom_point(data=plant.df2[(plant.df2$treatment=="burned"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
ggtitle("Burned") +
geom_line(aes(linetype= plant.type))+
ylab("Mass loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
compare.burned.mod.plot
#compare mass loss between sage and willow by fire treatment---unburned
m1.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
m2.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
m3.compare.unburned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
compare.unburned.AIC=AIC(m1.compare.unburned.gam, m2.compare.unburned.gam, m3.compare.unburned.gam)
compare.unburned.anova=anova(m1.compare.unburned.gam)
compare.unburned.anova
summary(m3.compare.unburned.gam)
gam.check(m3.compare.unburned.gam, rep=1000)
mean(plant.df2[plant.df2$plant.type=="sage",]$perc.loss)
mean(plant.df2[plant.df2$plant.type=="willow",]$perc.loss)
# model predictions
compare.burned.diff.mod<-plot_difference(
m1.compare.burned.gam,
series = plant.specific.mass..g,
difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod
#plot for the model output on rawdata
compare.unburned.mod.plot<-
plot_smooths(
model = m1.compare.unburned.gam,
series = plant.specific.mass..g,
comparison = plant.type
) +
geom_point(data=plant.df2[(plant.df2$treatment=="unburned"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
ggtitle("Unburned") +
ylab("Mass Loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
compare.unburned.mod.plot
##combine AIC
mod.comapre.decomp.rep<-rep(c("Burned: plant type with by-factor smooth",
"Burned: plant type with global smooth",
"Burned: global smooth",
"Unburned: plant type with by-factor smooth",
"Unburned: plant type with global smooth",
"Unburned: global smooth"))
mod.compare.decomp.column<- data.frame(mod.comapre.decomp.rep)
AIC.compare.decomp<-bind_rows(compare.burned.AIC, compare.unburned.AIC)
AIC.compare.decomp.mod<-as.data.frame(cbind(mod.compare.decomp.column, AIC.compare.decomp))
names(AIC.compare.decomp.mod)[1]="Model"
write.csv(AIC.compare.decomp.mod, "output/AIC models/AIC.comapre.decomp.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.AIC.pdf", width = 5.3, height = 2) # Export PDF
grid.table(AIC.compare.decomp.mod, rows=NULL)
dev.off()
compile raw plots and model-diff plots for final figures
decomp.compare.mod.plot<-plot_grid(
compare.burned.mod.plot+ theme(legend.position = "none"),
compare.unburned.mod.plot+ theme(legend.position = "none"),
ncol=2)
decomp.compare.mod.plot
ggsave("figures/decomp.comapre.mod.plot.pdf", height=4, width=12)Generate dataframe for model outputs
####Parametric output
decomp.compare.para.sums=as.data.frame(rbind(compare.burned.anova$pTerms.table, compare.unburned.anova$pTerms.table))
names(decomp.compare.para.sums)[1]="df/edf"
#create empty column for ref.df
decomp.compare.para.sums <- decomp.compare.para.sums %>%
# Creating an empty column:
add_column(Model= c("Burned: Plant type","Unburned: Plant type"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
decomp.compare.smooth.sums=as.data.frame(rbind(compare.burned.anova$s.table, compare.unburned.anova$s.table))
names(decomp.compare.smooth.sums)[1]="df/edf"
#create empty column for ref.df
decomp.compare.smooth.sums <- decomp.compare.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Burned: sage smooth","Burned: willow smooth","Unburned Model: sage smooth","Unburned Model: willow smooth"), .before="df/edf")
#merge the dataframes
decomp.compare.mod.allsums=rbind(decomp.compare.para.sums,decomp.compare.smooth.sums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.allsums.pdf", width = 6.2, height = 2) # Export PDF
grid.table(decomp.compare.mod.allsums, rows=NULL)
dev.off()
Comnine AIC and Mod sums for decomposition
library(gridExtra)
#merge all the dataframes
all.decomp.mod.sums=rbind(decomp.mod.allsums,decomp.compare.mod.allsums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.sums.pdf", width = 7, height = 3.7) # Export PDF
grid.table(all.decomp.mod.sums, rows=NULL)
dev.off()
#merge all the dataframes
all.decomp.mod.AIC=rbind(AIC.decomp.mod, AIC.compare.decomp.mod)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.AIC.pdf", width = 5.3, height = 3.7) # Export PDF
grid.table(all.decomp.mod.AIC, rows=NULL)
dev.off()
compile raw plots and model-diff plots for final figures
decomp.mod.plot<-plot_grid(
sage.mod.plot+ theme(legend.position = "none"),
willow.mod.plot+ theme(legend.position = "none"),
rel_widths = c(8,8), ncol=2)
decomp.mod.plot
ggsave("figures/decomp.mod.plot.pdf", height=4, width=12)compile raw plots and model-diff plots for final figures
all.decomp.mod.plot<-plot_grid(
decomp.mod.plot+ theme(legend.position = "none"),
decomp.compare.mod.plot+ theme(legend.position = "none"),
nrow=2)
all.decomp.mod.plot
ggsave("figures/all.decomp.mod.plot.pdf", height=8, width=8)